//ims_asc2tif.cpp

#include <string>
#include <vector>
#include <iostream>
#include <fstream>
#include "gdal_priv.h"
#include <sstream>
#include "wdatetime.h"
#include "filenameutil.h"
using namespace std;

#define MODE_0_4KM 0
#define MODE_1_24KM 1


bool readAscData(string infile, vector<unsigned char> &data )
{
	data.resize(6144*6144) ;
	std::ifstream infs(infile.c_str());
	std::string line;
    int numlines = 0;
    int datalines = 0 ;
	while (std::getline(infs, line))
	{
		++ numlines ;
		if( numlines >= 31 )
		{
			//datas
			if( line.length() == 6144 )
			{
				int iy = 6143 - datalines ;
				for(int ix = 0 ; ix < 6144 ; ++ ix )
				{

					int it = iy * 6144 + ix ;
					unsigned char val = 0 ;
					if( line[ix] == '1' ){
						val = 1 ;
					}else if( line[ix]=='2'){
						val = 2 ;
					}else if( line[ix]=='3'){
						val = 3 ;
					}else if( line[ix]=='4'){
						val = 4 ;
					}else{
						val = 0 ;
					}
					data[it] = val ;
				}
			}else
			{
				cout<<"data line size is not 6144 at line "<<numlines<<endl;
				return false ;
			}
			++ datalines ;
		}else if( numlines < 31 )
		{
			cout<<line<<endl ;
		}   
    }
    if( datalines==6144 )
    {
    	return true ;
    }else
    {
    	cout<<"date line count is not 6144, but "<<datalines<<endl ;
    	return false ;
    }
}

bool isUnpacked24km(string infile) {
	std::ifstream infs(infile.c_str());
	std::string line;
	std::getline(infs, line) ;
	infs.close();
	if( line.find("filehandle") != string::npos ) {
		return true ;
	}else{
		return false ;
	}
	
}

bool readAscData24Km(string infile, vector<unsigned char> &data )
{
	data.resize(1024*1024) ;
	std::ifstream infs(infile.c_str());
	std::string line;
    int numlines = 0;
    int datalines = 0 ;
	while (std::getline(infs, line))
	{
		++ numlines ;
		if( numlines >= 31 )
		{
			//datas
			if( line.length() == 1024 )
			{
				int iy = 1023 - datalines ;
				for(int ix = 0 ; ix < 1024 ; ++ ix )
				{

					int it = iy * 1024 + ix ;
					unsigned char val = 0 ;
					if( line[ix] == '1' ){
						val = 1 ;
					}else if( line[ix]=='2'){
						val = 2 ;
					}else if( line[ix]=='3'){
						val = 3 ;
					}else if( line[ix]=='4'){
						val = 4 ;
					}else{
						val = 0 ;
					}
					data[it] = val ;
				}
			}else
			{
				cout<<"data line size is not 1024 at line "<<numlines<<endl;
				return false ;
			}
			++ datalines ;
		}else if( numlines < 31 )
		{
			//cout<<line<<endl ;
		}   
    }
    if( datalines==1024 )
    {
    	return true ;
    }else
    {
    	cout<<"date line count is not 1024, but "<<datalines<<endl ;
    	return false ;
    }
}



bool readAscData24KmUnpack(string infile, vector<unsigned char> &data )
{
	data.resize(1024*1024) ;
	std::ifstream infs(infile.c_str());
	std::string line;
    int numlines = 0;
    int datalines = 0 ;
    int ix = 0 ;
    int iy = 1023 ;
	while (std::getline(infs, line))
	{
		++ numlines ;
		if( numlines >= 23 )
		{
			stringstream ss(line) ;
			int ival = 0 ;
			while( ss>>ival )
			{
				int it = iy * 1024 + ix ;
				if(ival==164){
					data[it] = 3 ;
				}else if( ival==165){
					data[it] = 4 ;
				}else{
					data[it] = ival ;
				}
				++ ix ;
				if( ix==1024 ){
					--iy ;
					ix = 0 ;
				}
			}
		}else if( numlines < 31 )
		{
			//cout<<line<<endl ;
		}   
    }
    return true ;
}

// PROJCS["unnamed",
//     GEOGCS["D_Unknown",
//         DATUM["unknown",
//             SPHEROID["unnamed",6378137,291.5053473491773]],
//         PRIMEM["Greenwich",0],
//         UNIT["degree",0.0174532925199433]],
//     PROJECTION["Polar_Stereographic"],
//     PARAMETER["latitude_of_origin",60],
//     PARAMETER["central_meridian",-80],
//     PARAMETER["scale_factor",1],
//     PARAMETER["false_easting",0],
//     PARAMETER["false_northing",0],
//     UNIT["metre",1,
//         AUTHORITY["EPSG","9001"]]]

const char PROJ_STR[]="PROJCS[\"unnamed\",GEOGCS[\"D_Unknown\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6378137,291.5053473491773]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",60],PARAMETER[\"central_meridian\",-80],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]" ;
const char PROJ_STR24km[]="PROJCS[\"unnamed\",GEOGCS[\"D_Unknown\",DATUM[\"unknown\",SPHEROID[\"unnamed\",6371200,291.5053473491773]],PRIMEM[\"Greenwich\",0],UNIT[\"degree\",0.0174532925199433]],PROJECTION[\"Polar_Stereographic\"],PARAMETER[\"latitude_of_origin\",60],PARAMETER[\"central_meridian\",-80],PARAMETER[\"scale_factor\",1],PARAMETER[\"false_easting\",0],PARAMETER[\"false_northing\",0],UNIT[\"metre\",1,AUTHORITY[\"EPSG\",\"9001\"]]]" ;

string convertname(string imsfullname,string outdir ){
	string bname = wFilenameUtil::baseName(imsfullname) ;
	int ymd7 =atof( bname.substr(3,7).c_str() ) ;
	int ymd8 = wDatetime::convertYmd7toYmd8(ymd7) ;
	char buff[32] ;
	sprintf( buff, "ims%d.tif" , ymd8 ) ;
	return outdir + "/"+string(buff) ;
}

int main(int argc, char* argv[] )
{
	cout<<"A program to convert IMS 4KM ASC file into geotiff file.2021-4-27"<<endl ;
	cout<<"usage: ims_asc2tif mode ims2001001.asc outdir "<<endl ;
	cout<<"input filename should be ims<yyyyddd>....asc format!"<<endl;
	cout<<"the asc file's data should be 6144x6144 , 31st line is data."<<endl ;
	cout<<"the 24km asc file's data should be 1024x1024 , 31st line is data."<<endl ;
	cout<<"no mode or mode=0 means 4km, mode=1 means 24km normal."<<endl ;
	cout<<"v1.0"<<endl ;
	cout<<"v2.0 support 24km asc data and unpacked data 164(3 ice) , 165(4 snow)"<<endl ;
	cout<<"v3.0 use outdir, add compress, auto check unpacked data."<<endl;

	GDALAllRegister() ;
	if( argc!=4 )
	{
		cout<<"params not 4."<<endl ;
		return 11 ;
	}

	int mode = 0 ;
	string infile , outdir  ;
	int outxsize = 0 ;
	int outysize = 0 ;
	double outreso = 0 ;

	mode = atof(argv[1]) ;
	infile = argv[2]; //"/home/hadoop/tempdata/IMS-ASC/ims1997035_24km_v1.1.asc" ;
	outdir= argv[3]; //"ims20210101.tif" ;	
	
	cout<<"input:"<<infile<<endl ;
	cout<<"outdir:"<<outdir<<endl ;
	cout<<"mode:"<<mode<<endl ;

	string outname = convertname( infile , outdir ) ;
	cout<<"outname:"<<outname<<endl;

	//第31行开始是数据
	vector<unsigned char> data;
	if( mode==MODE_0_4KM ){
		outxsize = 6144 ;
		outysize = 6144 ;
		outreso = 4000 ;
		bool isok = readAscData( infile , data) ;
		if( isok ==false )
		{
			cout<<"read asc data for 4km failed."<<endl ;
			return 12 ;
		}
	}
	else if(mode==MODE_1_24KM ) {
		//24km
		outxsize = 1024 ;
		outysize = 1024 ;
		outreso = 23684.997 ;

		bool unpack = isUnpacked24km(infile) ;
		if( unpack==true ) {
			cout<<"it is unpacked."<<endl;
			bool isok = readAscData24KmUnpack( infile , data) ;
			if( isok ==false )
			{
				cout<<"read asc data for 24km unpack failed."<<endl ;
				return 14 ;
			} 
		}else{
			bool isok = readAscData24Km( infile , data) ;
			if( isok ==false )
			{
				cout<<"read asc data for 24km failed."<<endl ;
				return 13 ;
			} 
		}
	}else{
		cout<<"unsupported mode."<<endl ;
		return 15 ;
	}

	{
		CPLStringList opt ;
		opt.SetNameValue("COMPRESS" , "Deflate") ;
		GDALDriver* driver = GetGDALDriverManager()->GetDriverByName("GTiff");
		GDALDataset* outds = driver->Create( outname.c_str(),outxsize, outysize,1,GDT_Byte, opt.List() ) ;
		if( outds ==0 )
		{
			cout<<"failed to create output file "<<outname<<endl ;
			return 13 ;
		}else
		{
			if( mode==MODE_0_4KM ){
				double trans[6] = {-12288000, outreso, 0 ,
							12288000, 0   ,-outreso } ;
				outds->SetGeoTransform(trans) ;
				outds->SetProjection(PROJ_STR);
			}else{
				double trans[6] = {-12126597, outreso, 0 ,
							12126840, 0   ,-outreso } ;
				outds->SetGeoTransform(trans) ;
				outds->SetProjection(PROJ_STR24km);
			}
			
			outds->GetRasterBand(1)->RasterIO(GF_Write,0,0, outxsize, outysize ,
				data.data(),outxsize,outysize,GDT_Byte,0,0,0) ;
			GDALClose(outds) ;
			cout<<"done"<<endl;
			return 0 ;
		}
	}

	return 0 ;
}